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(N ■ ABSTRACT 

I \ 

^ \ We present a 3-D NLTE Monte Carlo radiative transfer code that we use to 

^ \ study the temperature and ionization structure of Keplerian disks around Clas- 

sical Be stars. The method we employ is largely similiar to the Monte Carlo 
O \ transition probability method developed by Lucy. Here we present a simpli- 

\ fication of his method that avoids the use of the macro atom concept. Our 

\ investigations of the temperature structure of Be star disks show that the disk 

.!_H \ temperature behavior is a hybrid between the behavior of Young Stellar Object 

^ \ (YSO) disks and Hot Star winds. The optically thick inner parts of Be star disks 

\ have temperatures that are similar to YSO disks, while the optically thin outer 

parts are like stellar winds. Thus, the temperature at the disk midplane initially 
drops, reaching a minimum at 3-5 stellar radii, after which it rises back to the 
optically thin radiative equilibrium temperature at large distances. On the other 
hand, the optically thin upper layers of the disk are approximately isothermal 
— a behavior that is analogous to the hot upper layers of YSO disks. Interest- 
ingly, unlike the case of YSO disks, we find that disk flaring has little effect on 
the temperature structure of Be star disks. We also find that the disks are fully 
ionized, as expected, but there is an ionization minimum in the vicinity of the 
temperature minimum. The deficit of photoionization at this location makes it 
the most likely site for the low ionization state lines (e.g., Fe II) that produce 
the shell features observed in Be stars. Finally, we find that, despite the complex 
temperature structure, the infrared excess is well-approximated by an equivalent 
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isothermal disk model whose temperature is about 60% of the stellar tempera- 
ture. This is largely because, at long wavelengths, the effective photosphere of 
the disk is located in its isothermal regions. 

Subject headings: radiative transfer — stars: emission line. Be — circumstellar 
matter 



1. Introduction 

Circumstellar disks occur in a wide variety of astrophysical systems, ranging from Young 
Stellar Objects (YSOs) to rapidly rotating main sequence B stars (the Be stars) to highly 
evolved systems like cataclysmic variables, black holes, and perhaps asymptotic giant branch 
(ABG) stars and luminous blue variables (LBVs). In general, such disks tend to form 
whenever there is an accreting or outflowing circumstellar envelope in which the radial flow 
speed is small compared to the rotation speed. These disks reprocess the starlight, typically 
producing a strong infrared (IR) excess. Furthermore, these disks can be quite optically 
thick at some wavelengths, while being optically thin at others. Consequently, the radiation 
field at a given point can be quite non-local and anisotropic. In this paper, we develop 
a non-local thermodynamic equilibrium (NLTE) Monte Carlo radiation transfer code that 
may be applied to study differentially moving three-dimensional circumstellar envelopes. As 
a particular application of this code, we study the temperature structure of the disks around 
Be stars. 

The circumstellar environment of Be stars has been studied extensively in the past 
few decades (see Porter & Rivinius 2003, for a recent review). Both observational and 
theoretical evidence strongly suggests that most of the circumstellar matter is contained in a 
very dense, geometrically thin disk (Quirrenbach et al. 1997; Wood, Bjorkman, & Bjorkman 
1997). Since the stellar temperature is quite high, these disks are expected to be completely 
ionized and devoid of dust, so the disk opacity (primarily from hydrogen and free electrons) 
is much better understood than, for example, YSO disks (where the dominant opacity source 
is dust, whose size distribution, chemical composition, and even spatial distribution is rather 
uncertain). What is less well-determined is the detailed structure (density, temperature, 
ionization fraction, and velocity) of Be star disks. 

Kinematic information about the disk can be obtained from the disk emission lines. 
These line profiles arc often double-peaked with full widths (in the wings) that exceed the 
V sini of the star by about 20% (Hanuschik 1996), so the disk must rotate faster than the star. 
Based on the symmetry of the hydrogen Balmer line profiles and, especially the optically thin 
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Fell lines, it appears that the expansion velocity, Vr, of the disk must be much less than its 
rotation speed, (Hanuschik 1988, 1995, 1996, 2000). Studies that try to determine the ro- 
tation speed as a function of radius are less conclusive. Using a power-law parameterization, 
oc r^^ , Hummel & Vrancken (2000) find that they cannot distinguish between angular 
momentum-conserving (j = 1) and Keplerian (j = 1/2) velocities for optically thin lines, but 
for optically thick lines they find oc r~°'^^, favoring Keplerian orbital motion. Theoretical 
evidence for Keplerian rotation is provided by the observation of periodic variations in the 
shape of the double-peaked emission line profiles. During these variations, the amplitudes of 
the violet and red emission peaks change in an antisymmetric fashion (so-called V/ R varia- 
tions). Such behavior can be explained theoretically by a precessing one-armed density wave 
(Okazaki 1991), but on theoretical grounds such density waves require a Keplerian velocity 
structure (Okazaki 1991). We conclude, therefore, that the disk has a Keplerian orbital 
velocity with v^. <^v^. 

The IR excess of Be stars arises from free-free (f-f) emission in the disk (Gchrz. Hackwell, 
& Jones 1974). Since the f-f opacity increases with wavelength, the slope of the IR excess can 
be used to probe the radial density structure of the circumstellar envelope (Wright & Barlow 
1975; Lamers & Waters 1984). Assuming an isothermal disk with a radial power law for 
the density, p oc r~", this technique gives density exponents in the range n = 2-3.5 (Waters 
1986; Cote Sz Waters 1987; Waters, Cote, Sz Lamers 1987). It should be noted however, 
that this result depends on the radial structure of the disk temperature (see CassineUi & 
Harmann 1977, which shows how the slope of the IR excess depends on the temperature 
gradient) . 

Although quite sophisticated radiative equilibrium models have been developed that 
can determine the 2-D temperature structure of YSO disks and envelopes (e.g., D'Alessio et 
al. 1998; Wood et al. 2002; DuUemond 2002; Whitney et al. 2003; Walker et al. 2004), we 
are only now beginning to understand the temperature structure of Be star disks. Models 
of spherically expanding stellar winds suggested that, as a consequence of the Cayrel effect 
(Cayrel 1963), the wind ionization fractions would "freeze" (i.e., be constant) at large radii 
and that the temperature would become isothermal with a temperature of about 0.8Teff 
(Klein & Castor 1978, however later work by Drew 1989, who included cooling by heavy 
elements, shows the temperature is not particularly isothermal, slowly faUing to less than 
60% TeG at large distances). In any case, the usual assumption for the disk temperature was 
to hope for a similar situation and to assume that the disk was isothermal with a tempera- 
ture Ta — 0.8Teff (e.g.. Waters 1986). The first attempts to calculate the disk temperatures 
were performed by Millar & Marlborough (1998), who implemented the radiative equilibrium 
condition within the Poeckcrt & Marlborough (1978) code. Note however, that the Poeck- 
ert & Marlborough code does not perform a full two-dimensional solution to the radiative 
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transfer problem; it instead employs several approximations to estimate the hydrogen level 
populations. In subsequent papers, Millar & Marlborough (1999a) approximated the effects 
of diffuse emission by employing the "on-the-spot approximation", while Jones, Sigut, & 
Marlborough (2004) included cooling by Fe line emission. In general, they find that the disk 
temperature in the mid-plane is quite low, while the upper layers of the disk are hot. This 
fundamental work shows that the disk is certainly not isothermal, bringing into question the 
Waters et al. density results. For this reason, Millar & Marlborough (1999b) calculated the 
temperature for the Waters (1986) density model. They concluded that the average disk 
temperature is a factor of 2-3 lower than that assumed by Waters. An additional effect (not 
addressed by Millar & Marlborough 1999b) is that the non-isothermal temperature structure 
may alter the slope of the IR continuum, thereby changing the estimates of the radial density 
exponent, n. 

In this paper, we build upon and extend the results of Millar & Marlborough (1998, 
1999a,b) to study more fully the temperature and ionization structure of the disk. In partic- 
ular, we investigate what physically determines the temperature behavior of the disk, and we 
compare this behavior to that of YSO accretion disks and stellar winds. We also examine the 
effects of a non-isothermal disk temperature on the IR spectral energy distribution (SED) of 
Be stars. To do so, we have developed a Monte Carlo radiative transfer code (presented in 
section 2) that fully solves the 3-D NLTE problem and radiative equilibrium temperature. 
In section 3, we present test calculations used to verify both the method and the computer 
code. The disk model used in this paper is described in Section 4. In Section 5, we investigate 
the dependence of the temperature, hydrogen level populations, and ionization fractions on 
model parameters such as the disk density and geometry. Similarly, we investigate how the 
temperature structure alters the slope of the IR excess. Finally, we summarize our results 
in Section 6. 



2. The Monte Carlo Code 

One of the greatest assets of the Monte Carlo (MC) method is its flexibility and ease of 
implementation, especially as the complexity of the system increases. Although traditional 
radiative transfer methods are much more computationally efficient for simple configurations 
(e.g. plane-parallel or spherical geometries), for more complex situations (e.g., 3-D geome- 
tries, anisotropic scattering, and complex frequency redistribution processes) the MC method 
is frequently the only viable alternative. Furthermore, the relative efficiency of Monte Carlo 
methods increases as the number of spatial dimensions, frequencies, and anisotropy of the 
radiation field increases. 
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Recently, Lucy (2002, 2003, 2005, hereafter L02, L03, and LOS) developed a Monte 
Carlo transition probability method to solve the statistical equilibrium problem. His method 
employs indivisible energy packets to describe the exchange of energy between photons, 
macro atoms (an ensemble representation of the atomic excitation throughout a finite volume 
element), and free electrons. The energy packets contain either radiative energy (r- packets) 
or kinetic energy (/c-packets). When an energy packet (either radiative or kinetic) is absorbed 
by a macro atom, the macro atom randomly performs internal transitions (using an effective 
transition probability array) between its energy states until it emits a new energy packet 
(either kinetic or radiative); see Figure 1. Since any photon absorption is eventually followed 
by emission of an equal energy photon packet, this method automatically enforces radiative 
equilibrium (in the rest frame of the gas). The key ingredient of this method is the effective 
transition probability array, and L02 shows how to choose these probabilities so that the 
internal transitions of the macro atoms will reproduce the statistical equilibrium equations. 
Since the transition probabilities (effectively they are the net rates between levels) depend 
on both the values of the mean intensity in each line and the level populations, their solution 
requires iteration. 

L02 shows that this technique correctly reproduces the line and continuum emissivities, 
the "A-iteration" converges, and on average there are only about two internal transitions 
per photon- macro atom interaction. In subsequent papers (LOS and LOS), Lucy applied 
this method to the non-trivial problem of calculating the hydrogen level populations and 
electron temperature in the expanding shell of a type II supernova. He also studied the 
convergence behavior and demonstrated that the correct solution can be achieved with a 
reasonable amount of computing time. 

In 2002-2004, we independently developed a method for solving both the statistical 
and radiative equilibrium in circumstellar winds and disks of hot stars. Our method is 
largely similar to the one described in L03, in that we measure the rates (radiative, heating, 
and cooling) using essentially identical methods — path length sampling, (see Lucy 1999, 
hereafter L99), so rather than reproducing those formulae here, we simply refer the reader 
to sections 5 and 6 of LOS and section 2 of LOS. However our method differs from Lucy's in 
one key area: the photon-atom interactions. In the following, we present a short description 
of our method, focusing on the differences with Lucy's method. 

The essence of any method for solving the radiative transfer problem is to correctly de- 
scribe the radiation field everywhere in the so-called interaction region (the physical medium 
wherein the radiation propagates). In the MC method, this is accomplished statistically by 
simulating the random propagation of N photon packets (PPs) through the medium, where 
they are scattered, absorbed, and reemitted until they eventually leave the interaction re- 
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gion. A convenient choice is to emit equal-energy monochromatic PPs from the star (as a 
consequence, packets of differing frequencies contain different numbers of physical photons). 
The condition of radiative equilibrium then implies that whenever a PP is absorbed, it must 
be reemitted at the same location with the same energy (in the rest frame of the fluid; in 
the stellar frame, the packet energy changes, owing to the doppler shift between the different 
reference frames). This device automatically ensures flux conservation; i.e., the flux is ev- 
erywhere divergence free (Lucy 1999; Bjorkman & Wood 2001). Although the (rest frame) 
packet energy does not change, the new direction and frequency of the PP must be chosen 
to reproduce the emissivity ji, of the medium at that location. 

To solve the statistical equilibrium rate equations, the following procedure is employed. 
Initially, the interaction region is divided into a number of cells. In each cell, the gas state 
variables (density, velocity gradient, kinetic temperature, and level populations) are assumed 
to be constant. The simulation starts with an initial estimate of the state variables for each 
cell. While the simulation is performed, the radiative rates, as well as the heating rate, of 
the gas are sampled as photons pass through each cell. At the end of the MC simulation, 
the accumulated radiative and heating rates are used to update the cell level populations 
and kinetic temperature. A new simulation then is performed using the updated values of 
the state variables. This procedure is repeated until the state variables converge to their 
equilibrium values for all cells. 

From the above description, it is evident that two key points arc necessary for the method 
to succeed. First, there must be an efficient way to sample all the radiative transitions and 
heating rates. An elegant and efficient solution for this was proposed by L99 and further 
developed in L03. In section 2.1.2, we describe how to more efficiently sample the bound- 
bound (b-b) rates. The second critical point is that the MC simulated radiation field has to 
statistically reproduce the "true" radiation field of the system. Because every PP absorption 
is immediately followed by in situ reemission, this condition reduces to sampling correctly 
the frequency and direction of the reemitted PPs. 

Lucy's method (fig. 1) simulates the interaction of the radiation with matter by tracking 
transitions of his macro-atoms. When a PP (r-packet) undergoes an absorption event (de- 
termined by the opacities), it will either activate a macro atom (for instance, through b-b or 
bound-free transitions) or be converted into a fc-packet (for example, via f-f absorption). A 
series of internal transitions (controlled by the MC transition probabilities, L02) eventually 
leads to the conversion of the energy back into an r-packet. Thus, Lucy's method is a scheme 
for sampling the frequency of the reemitted PPs. 

The approach we adopt in this paper differs from Lucy's in that we do not use his MC 
transition probabilities to sample the frequency of the reemitted PP. Instead, we bypass 
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his macro-atoms and directly sample the frequency from the total emissivity of all possible 
emission processes. Thus our emission process is completely unlinked from the absorption 
process. For example, a photon may be absorbed by a b-b transition and then reemitted by a 
f-f process. This may seem peculiar, but all that matters is that at the end of the simulation, 
we produce the correct total luminosity and frequency distribution for each process. 

For a given cell j, the total luminosity of a given process k is 

Lj= / dV / Airf^diy. (1) 
JVj Jo 

where is the emissivity of process k and the first integral is performed over the cell volume 
Vj. When a PP is absorbed in cell j, we select the physical process k that reemits the PP, 
determined by the minimum value of k satisfying 

k 

e<^L5/L„ (2) 

i=l 

where ^ is a uniform random number in the interval (0, 1), and the total luminosity of the 
cell 

^.^E^5- (3) 

i 

Once the physical process is chosen, the frequency of the emitted PP is sampled using 
standard MC techniques (see section 4 of L03). 

From the point of view of computer coding, our method has the advantage of being 
easier to implement (but one should note that the implementation of macro-atoms does not 
require a particularly large coding effort). Since we dispense with the macro-atom device, we 
do not need to track their internal transitions, which in principle increases the relative speed 
of our code (but, this may not be a large effect since there were only two internal transitions 
per interaction for the case studied by L02). Similarly, there is no need to convert energy 
packets back and forth between r-packets and A;-packets. On the other hand, it is not 
clear which method will converge faster. Both methods use the previous line intensities 
to calculate their reemission probabilities, and the two methods are rigorously equivalent 
for a system already in equilibrium. However, out of equilibrium, the two methods differ. 
Our emissivities are constant during a given iteration, but Lucy's method uses conditional 
transition probabilities, which may somewhat adapt the emissivities to the new values of the 
line intensities. We suspect, therefore, that Lucy's method might have a somewhat better 
convergence rate. 

Finally, we point out that even using the efficient intensity estimators of L03 and LOS, the 
computational demand for 2-D and 3-D problems is enormous, requiring a parallel-computing 
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cluster to make the run times practical. Fortunately, the structure of the algorithm is well- 
suited for parallel computation. The problem is "embarassingly parallel" — one simply 
divides the total number of photon packets into batches that run independently on each 
processor. For the 2-D models shown in Sections 5, we used a 96-node Beowulf cluster (with 
1.7 GHz Athlon CPUs). One iteration for the level populations and kinetic temperature 
requires about 15 minutes on this cluster (using N = 3 x 10^ PPs). Since the number of 
iterations required is about 10-15, one model takes a few hours to complete. 



A common problem for full-spcctrum MC codes is that the radiation field is typically 
under-sampled at frequencies where the mean intensity is low. This happens because the 
number of PPs per frequency interval is proportional to the mean intensity. As a result, 
calculating the emergent flux at those frequencies requires a very large number of PPs for 
the simulation. Similarly, calculating the emergent flux for a very narrow frequency range 
(as is the case for line profiles) can be problematic because the number of PPs in that range 
is small. 

A related issue is determining the b-b rates required to solve the statistical equihb- 
rium equations, because sampling each transition requires a large number of PPs in a small 
frequency range. For instance, suppose that the mean intensity at a given point of the in- 
teraction region is given by a Planck function. The ratio of resonant Ha PPs to the total 
number of PPs is approximately 



where Wmax is the maximum macroscopic velocity of the gas and a is the Stcphan-Boltzmann 
constant. For fmax = 1000 kms^^ and T = 16000 K, this ratio is about 10^^, which 
illustrates the low number of PPs available for sampling the b-b rates. If one considers that 
the main source of exciting PPs comes from the external radiation source (the star), and 
that those PPs would be spht over the very large number of cells necessary for 2-D or 3-D 
problems, it is easy to see that a huge number of PPs would be necessary to correctly sample 
the b-b rates. The solution to this problem is simple: limit the frequency range of the MC 
simulation to the vicinity of the line of interest, and sample each b-b transition separately. 

Consider a star surrounded by an envelope that is divided into M cells. The total energy 
emitted by the star per unit time in the frequency range uq and i/i is 



2.1. Simulations for a Small Frequency Range 



N 



(4) 




(5) 
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where is the outward flux at the steUar surface. Similarly, the luminosity of process k 
in cell J, L*^, can be calculated from eq. (1), using the appropriate frequency range in the 
integral. The total energy emitted by the cell, Lj, is calculated from eq. (3), and the total 
energy emitted by the envelope is given by 

M 

L^nv = ^Lj. (6) 

J=l 

Given the energies calculated above, the MC calculation proceeds step-by-step as follows: 

1. A new PP is generated. If ^ < L-^/ (L^ -|- Lenv) the PP is emitted by the star (go to 
step 2), and if C > -^*/ (-^* + -^env) the PP is emitted by the envelope (go to step 3). 

2. The position on the stellar surface is sampled randomly, the PP frequency is chosen 
from the stellar flux, JF+, over the relevant frequency range, and the propagation 
direction is sampled from the limb darkening law. Proceed to step 6. 

3. The cell j that emits the PP is sampled using the relation ^ < J2i=i Li/^env 

4. The process k that emits the PP in cell j is sampled using eq. (2). 

5. The position of the PP inside the cell is chosen randomly, and the PP propagation 
direction is chosen accordingly to the process of emission. For example, the direction 
distribution for line emission is given by the Sobolev escape probability. 

6. The PP is propagated through the medium using standard MC techniques until it 
leaves the envelope. If the PP is absorbed or hits the star, it is destroyed. Return to 
step 1 until all N PPs are generated. 



2.1.1. Detailed Spectra 

The above procedure can be used for the calculation of emergent spectra, both where 
the flux is low (e.g., the IR flux in hot stars) and when the frequency range is small (e.g., 
a detailed calculation of the shape the balmer jump). Recall that in the radiative equi- 
librium scheme, a PP may be absorbed and rccmitcd many times before escaping — a 
time-consuming process — while in the above procedure the PP is simply destroyed. Thus 
in addition to enabling the calculation of multi-frequency radiative transfer in arbitrary 
frequency ranges, the above procedure also has the advantage of performing much quicker 
spectral syntheses. Furthermore, it automatically includes anisotropic scattering and polar- 
ization effects that are difficult to include in traditional post-processing routines. 
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2.1.2. Sampling Bound- Bound Rates 

Finally, wc use the above technique to sample the b-b rates in the following way. 
We first run a standard radiative equilibrium simulation to calculate all the radiative and 
heating rates, apart from from the b-b rates. Then a series of simulations are performed 
to sample individually the rates for each b-b transition. By choosing a frequency range 
[^'l(1 — 2vniax/c), 1^(1 -|- 2i>inax/c)], where i/L is the frequency of the transition (i.e., centered 
on the line of interest), we ensure that we properly account for the coupling to all possible 
processes, including other emission lines. 

3. Numerical Tests 

Here we verify our code results and test the convergence of our algorithm. We present 
two test cases: 1) a blackbody cavity, which has a simple analytic solution, and 2) the 
non-trivial problem of a spherically symmetric expanding stellar wind. 

3.1. Blackbody Cavity 

This test consists of a spherical cavity whose boundary emits a blackbody spectrum of 
a given temperature. The cavity is filled with hydrogen with constant density throughout. 
Each PP emitted by the cavity wall undergoes a random series of scatterings and absorptions 
followed by in situ reemissions as described in section 2, but eventually the PP is absorbed 
(and destroyed) by the blackbody cavity wall, whereupon another PP is emitted until the 
total number PPs have been run. 

A blackbody cavity is an especially convenient test because the gas should come to 
complete thermodynamic equilibrium with the wall. Thus the state variables of the gas 
depend solely on the temperature of the cavity's wall, Tcav, and the level populations and 
ionization fraction should be given by the Boltzmann-Saha distribution. Furthermore, the 
local radiation field should be a Planck function, and all radiative rates should attain their 
LTE values. Finally, the kinetic temperature of the electrons should also be Tcav 

We used this setup to test several aspects of our method. A trivial but important 
situation is to set the initial level populations and electron temperature to their correct 
equilibrium values; i.e., Tg = Tcav and rii = n*, where the asterisk indicates LTE level 
populations. In this case, one expects that, if the rates are being properly sampled, the 
state variables should oscillate around the LTE values as a result of MC statistical sampling 
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errors. Furthermore, the magnitude of the oscillation should decrease approximately^ as 
1/\/N. This test is also important because it demonstrates the stability of a converged 
solution. The results for 30 iterations with initial conditions T^av = 20000 K and hydrogen 
number density Uh = 10^^ cm~^ are shown in Figure 2. The simulation in the left panel 
employed = 10^ PPs, while the right panel had N = 10^. As expected, the results for Tg 
and rii oscillate around the correct equilibrium values, and the magnitude of the oscillation 
decreases as N increases. The R.M.S. relative error 

1/2 

(7) 

where S is the state variable of interest, is indicated for hydrogen levels 1 and 10, as well 
as the electron temperature. Note that the R.M.S. errors decrease roughly by the expected 
amount (VlO^). 

In the next series of tests, we varied the initial values of the state variables to verify that 
the results did converge to their correct values. As an example. Figure 3 shows 30 iterations 
of the ground level population for three values of N (10^, 10^, and 10^). The population 
was increased initially by a factor of 10^, and it converges to the correct value (shown as a 
straight line in the figure) after about 10 iterations. More importantly, the level population 
converges even in the presence of the very large sampling errors for the A^ = 10^ case. 

The high precision we were able to obtain for the level populations and electron tem- 
perature in the blackbody cavity test verifies the validity and the correctness of our imple- 
mentation of the method described in the last section. In addition, it shows that, as already 
demonstrated by L02, the random errors intrinsic to the Monte Carlo method do not prevent 
the convergence of the state variables. It is important to emphasize this last remark because, 
when the method is applied to 2-D or 3-D geometries, one may easily run into situations for 
which given parts of the simulation may be under-sampled. An example of such a situation 
is the equatorial region of optically thick disks, which is shielded from the stellar ionizing 
radiation by the disk's upper layers. 

Finally, we would like to comment that the blackbody cavity test has proven extremely 
useful during the development phase of the code. In particular, detailed balance occurs in 
LTE, so one can turn on and off individual processes (e.g., f-f emission and absorption), 
and the results should not change. This permits one to test separately each part of the 
implementation. Similarly, sampling the local radiation field and looking for departures 



riit 



H "it 



S*f/{S*) 



^ Since there is a nonlinear relationship between the rates and the state variables, the statistical error will 
not follow exactly a Poisson distribution (unless the fluctuations are small). 
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from a Planck function is especially useful for diagnosing errors. 



3.2. NLTE Test 

To assess the NLTE results for our method and to verify the Sobolev escape probabilities 
for the spectral lines (which modify the rate equations), we applied our code to the 1-D 
problem of a spherical expanding wind around a hot star, and compared our results to those 
given by the well-estabhshed code of MacFarlane et al. (1993). We were able to match 
MacFarlane's results for a wide range of stellar temperatures, wind density, and velocity. An 
example comparison is given by Figure 4, which shows the n = 2 hydrogen level populations 
for a B2V star (Tejf = 20000K, = 6.6Rq), surrounded by a stellar wind with a /9-law 
velocity 

Vr ^Voo(^l- , (8) 

with Voo = 1600 km s"^ and (3 = 0.8. For all mass-loss rates shown (logM = —8.7, —7.7 
and —6.7), the difference between the results is typically less than 0.1%. 



4. Disk Model 

As discussed in the introduction, observational evidence indicates that the circumstellar 
disks of Be stars are geometrically thin, optically thick disks, whose outflow velocities are 
small in comparison to their rotation speeds. Dynamically this requires that the disk be pres- 
sure supported in the vertical direction and centrifugally supported in the radial direction. 
For this reason the disk structure is most likely to be that of a Keplerian hydrostatically 
supported disk. How material is injected at Keplerian speed into this disk is still unclear, but 
after it is injected, the gas will spread outward (and inward if it is injected at intermediate 
radii), owing to viscous effects (Lynden-Bell & Pringle 1974). 

Determining the self-consistent structure for a Keplerian disk is a somewhat complicated 
problem, because the disk temperature controls its geometry (via the hydrostatic equilibrium 
equation), which in turn determines the heating and hence temperature of the disk (see 
Kenyon & Hartmann 1987). Rather than solving the hydrostatic equilibrium equation (which 
we plan to investigate in a subsequent paper), we instead assume that the disk geometry is 
a known function. For our present purposes — an initial investigation of the temperature 
and ionization structure of the disk — we feel this is a reasonable approach. 

Studies of pre-main-sequence accretion disks (see Hartmann 1998, for a review) show 
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that the disk pressure scale height H — {a/v^)zu is small (since the orbital speed, v^, is 
typically much greater than the isothermal sound speed a), so the disk is geometrically thin, 
and the disk temperature typically falls (roughly) as a power-law with radius. Assuming the 
disk is vertically isothermal, the disk density structure is given, in cylindrical coordinates 

(OT,^,0),by ^ ^ 

p(^,.)=Po(§) exp(-^) (9) 

(see, for example, Bjorkman 1997). For steady-state isothermal outflow, the radial density 
exponent n — 3.5 (Lee et al. 1991; Porter 1999). The vertical structure is the familiar 
Gaussian result for pressure-supported disks with a scaleheight 

H{w)^Hoi^^J , (10) 

where 

Ho = —R., (11) 

and T;crit = \^GM/E^ is the critical velocity of the star. The disk flaring exponent (3 
depends on how the disk temperature falls with radius. For optically thick, infinitessimally 
thin reprocessing disks, T oc w^^^^, which gives P = 9/8 (Kenyon & Hartmann 1987, but 
see D'Alessio et al. 1998), while for optically thin stellar winds, the temperature is roughly 
isothermal (Klein & Castor 1978), so /3 = 1.5. 

The disk flaring parameter /3 controls the shape of the disk surface and hence its heating 
versus radius. As such, it is the primary parameter, along with density, that controls the 
temperature and ionization structure of the disk. Consequently for this paper, we adopt the 
above structure for the circumstellar disk, and let f3 and po be free parameters. Some models 
of Be star disks assume the disk is truncated at some point, perhaps by a binary companion, 
so we also assume the disk has a flnite outer radius i?d, which we set to lOOOi?*. 

To complete our disk model, we have to specify the hydrogen model. We adopt a 25-level 
model of the H atom, where the first 12 levels are explicit NLTE levels, and the remaining 
13 upper levels are implicit levels with LTE populations (with respect to the continuum). 
Each level has a consolidated statistical weight Qn = Iv?'. 

For the description of the central star we have to specify its efi^ective temperature and 
surface gravity, which defines the Kurucz model atmosphere (Kurucz 1979, 1994) that we 
use as our input spectrum. We have to specify, in addition, the stellar radius and mass, 
which determine the critical velocity and scaleheight Hq. Previously, we modeled the disk of 
C Tauri, B3IV (Wood et al. 1997), so we arbitrarily chose this spectral type for our stellar 
parameters (see Table 1). 
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4.1. Grid Definition 

As described in section 2, the MC calculation is performed for a discretized model of the 
disk with each of the model unknowns (rates, temperature, level populations) being constant 
throughout a given cell. Thus, it is essential that the grid be constructed to minimize (as 
far as practical) the variation of these unknowns within the cell. 

The disk is divided into radial shells, each of which are divided into cells in latitude 
and, if needed, into cells in azimuth. We found a total of 40 cells in radius and 40 cells 
in // to be a good compromise between usage of computer memory, speed of the code, and 
accuracy of the result. The radial spacing is defined so that each cell has the same radial 
electron optical depth. This approach has the drawback that the last cell becomes very large, 
so we further subdivide the last cell to properly sample the outer parts of the disk. The 
//-spacing is defined in such a way that each cell represents an equal linear step in density, 
and, consequently, the cells are smaller close to the midplane and become larger for larger 
z. Also, the cells become larger with increasing radius because of the disk flaring. The disk 
is truncated in the 2;-direction when the density becomes 10~^ smaller than the value at the 
midplane. 

5. Results 

In this section we explore the effects of the disk free parameters (po and (3) on the 
electron temperature and hydrogen level populations. We also investigate the effects of the 
temperature distribution on the IR excess. The parameters for all the models shown in the 
subsequent sections are listed in Table 2. 

5.1. Ionization and Level Populations 

The n = 1 and 2 level populations of hydrogen for models 01 and 04 (the density ex- 
tremes) are shown in Figure 5, as a function of radius. Each panel shows the level populations 
at different latitudes as indicated (0° corresponds to the disk midplane). The approximations 
used in Millar & Marlborough (1998, 1999a) make it difficult to perform a direct comparison 
between our results and theirs; however, they are qualitatively similar. 

In general, the disk is almost completely ionized, with the exception of the midplane. 
However, in none of the models does the midplane become completely neutral. For model 04 
(the highest density disk), the neutral hydrogen fraction at the midplane has two maxima. 
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one at fa GR^, for which ni — 20% and the second ed, zu ^ 50R^ where ni — 44%. For 
model 01 (the lowest density model), the two maxima are also present, but closer to the star 
(rii = 1% at zu ^ 3i?* and ni = 6% at w ^ 40R^). These minima in the hydrogen ionization 
fraction indicate the location where the photoionizing mean intensity is smallest. The first 
minimum is also near the temperature minimum of the disk (see section 5.2), so this position 
is likely to be the location of the low ionization state lines (like Fe II) that form the shell 
features in Be spectra. 

Another interesting feature of Figure 5 is the crossing of the low latitude n = 1 curves. 
At small radii, the lower density model (01) becomes more neutral than the higher density 
model (04). This apparent incongruity is a result of the higher density model being much 
hotter near the star than the lower density model. 

5.2. Temperature 

The temperature structure of model 02 is shown in Figure 6. The leftmost panel is the 
temperature in the vicinity of the star {1 < to/R^ < A), the center panel is the temperature 
for intermediate distances (1 < zu/R-/, < 10), and the rightmost panel is for large distances {w 
up to lOORi,). In this figure, which generally is representative of all models, we observe two 
regions with very distinct temperature behaviors: the midplane region (between the dashed 
curves, defined by \z\ < H), and the upper/lower layers of the disk {\z\ > H). In the upper 
layers, the temperature is nearly constant, averaging 10100 ± 1300 K (about 50% of Teff). 
This region is optically thin for shorter wavelengths (typically A < 10 //.m), so it contributes 
very little to the disk emission at these wavelengths; however, it becomes optically thick for 
longer wavelengths. Thus, this layer controls the long wavelength IR excess. 

For the midplane region, the temperature structure is much more complex. Beginning 
at the base of the disk, the temperature is high {k. 15000 K), quickly drops to a minimum 
temperature of 7000 ¥^ at w 4i?^, and rises back, reaching ^ 11000 ¥^ at w lOi?^, 
beyond which it becomes roughly isothermal. Other models show similar behavior, but the 
depth and location of the temperature minimum depends on the disk density (or equivalently, 
optical depth) . Figure 7 shows the mid- plane temperature as a function of increasing density 
(models 01-04). Note that as the density increases, the temperature minimum becomes 
deeper and occurs at larger radii. 

Since the disk is strongly non-isothermal (unlike a stellar wind), it is important to 
understand what determines its temperature structure. In accretion disks around young 
stellar objects (YSOs), the temperature typically falls as a power law with radius, T oc r~"*. 
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with an exponent m in the range 1/2-3/4, depending on the disk flaring. This is because 
the surface of a very optically thick disk acts very much like a blackbody that reprocesses 
the incident starhght. Since the star iUuminates the surface at an increasingly obhque angle 
with radius, the temperature falls faster than the usual r^^/^ geometrical dilution of the star. 
For an infinitessimally thin, flat (i.e., not flared) blackbody disk, the temperature falls as 
r~^/^, but close to the star a detailed integration gives 



where is the temperature of the radiation that illuminates the disk (Adams, Lada, & Shu 
1987, hereafter, ALS). 

Using this function, we flt the mid-plane temperatures (close to the star) for each model. 
The resulting curves for models 01-04 are shown in Figure 7 (thick lines) , while the values of 
along with vndep (the outer radius of the fitting interval) are listed in Table 2. We see that 
eq. (12) reproduces the initial temperature decrease quite well. However, is larger than 
Teff for all models (except the lowest density model 01), and the larger the density, the larger 
the corresponding T^. It is evident, therefore, that this behavior is due to back-warming of 
the star (photons that hit the star are reemitted locally to simulate this effect). This result 
reveals a small inconsistency in our model: the back-warming is sufficiently large that the 
stellar photospheric temperature should be increased near the equator. We did not attempt 
to fully model this effect; we only reemitted the photons without raising the temperature. 

Another interesting feature of Figure 7 is that, as the density increases, the departure 
point, ci7dcp (the location where the MC temperature curve departs from eq. [12]) shifts to 
larger radii. This suggests that Wdep is controlled by the disk optical depth. Indeed, we 
found that zu^ep correlates very well with the vertical electron scattering optical depth 



where ut is the electron scattering cross section, 7 is the number of electrons per ion, jj, is 
the mean molecular weight, and mn the mass of the hydrogen atom. The values of Te(rodep) 
for each model are fisted in Table 2. We conclude that the temperature stops falfing because 
the disk becomes (vertically) optically thin. Once the disk becomes optically thin, ionizing 
photons from high stellar latitudes begin to penetrate the disk, heating it back up to the 
optically thin radiative equilibrium temperature at large radii. Beyond this location the disk 
becomes roughly isothermal, behaving much like a optically thin stellar wind. 

The behavior of the midplane temperature can be summarized as follows: close to the 
star, the temperature is described by a flat blackbody reprocessing disk, eq. (12), with 




(12) 




(13) 
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Tes for low po and T^, > Tes for larger po- When the vertical electron optical depth 
drops below Te ~ 0.075, the temperature curve departs from eq. (12) and quickly returns to 
the optically thin radiative equilibrium temperature (about 60% T^g) for w > 10-30i?^. 

We end this section with a comparison between our results and those of Millar & Marl- 
borough (1999a,b). For the optically thin parts of the disk, our results agree — we also find a 
nearly constant temperature of about 60% Tcfj. For the midplane, however, our results differ. 
Their temperature in the outer parts of the disk is in general lower than ours, while in the 
inner parts, their results are often higher. In addition, their mid-plane temperature profiles 
initially drop much more quickly than ours, and they sometimes have a local maximum at 
intermediate radii that is not seen in our results. Nonetheless, we do agree qualitatively that 
there is temperature minimum at small radii in the disk mid-plane. 

5.3. Disk Flaring 

Let us now investigate the effects of the disk flaring parameter (3 on the temperature 
structure. Experience with the flared disks of YSOs shows that the flaring controls the 
amount of heating versus radius, which determines the temperature versus radius (e.g., 
Kenyon & Hartmann 1987). In general, one finds that the smaller the fiaring, the steeper 
the temperature gradient. 

We do not, however, observe this in our results. Figure 8 shows the temperature as 
a function of radius for models 05, 06, 07, 08, 09 and 02, which have (3 ranging from 1.0 
(wedge-shaped disk) to 1.5 (radially isothermal disk). The left side shows the results for 
three different latitudes. We see that the mid-plane temperature is very similar for all the 
models, but for the other latitudes, the temperatures are quite different. This difference, 
however, is a result of the different scalcheights probed by each latitude ray. If instead, we 
plot the temperatures at fixed fractions of the scaleheight above the mid-plane (right side), 
the curves are all very similar. We conclude that the disk fiaring does not infiuence the disk 
heating versus radius. 

There are probably two reasons for this: 1) at small radii, the disk surface is much less 
obliquely illuminated, and the direction perpendicular to the disk surface is insensitive to f3 
(this is why the ALS model works well for this region); and 2) the outer parts of the disk 
are optically thin, so they heat up to the optically thin radiative equilibrium temperature, 
regardless of the disk fiaring. 
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5.4. 



IR excess 



The spectral energy distribution for Model 02 in the 0.08-1000 /xm wavelength range 
is shown in Figure 9. Interestingly, in addition to the usual IR excess, there is a small UV 
excess (in the Balmer continuum) for pole-on viewing angles. Just longward of the Lyman 
jump, this UV excess is a result of electron scattering in the disk, which redirects some of 
the stellar flux into the polar direction. As the wavelength increases, bound-free emission 
begins to dominate, so the pole-on excess increases until the Balmer jump occurs. Beyond 
the Balmer jump, the f-f emission (and opacity) becomes increasingly important. As the 
envelope becomes optically thick, the disk develops a pseudo-photosphere that grows in size 
with increasing wavelength, producing the large IR excess. Since the disk is geometrically 
thin, the area of the pseudo-photosphere is smaller viewed edge-on than pole-on. As a result, 
the edge-on IR flux is about 10-50% of the pole-on flux, depending on wavelength. Note 
that the slope of the IR excess also depends on inchnation angle. 

As mentioned in the introduction, the slope of the IR excess has been used by Waters 
(1986) to observationally measure the radial density exponent for the circumstellar envelopes 
of Be stars. However, as demonstrated in section 5.2 (also see Millar & Marlborough 1998), 
the temperature varies throughout the disk by more than a factor of two, which raises 
questions about the validity and accuracy of Waters' results. 

To isolate and assess the effects of the temperature structure, we compare the IR excess. 



of the radiative equilibrium model to that of an equivalent isothermal model (i.e., one with 
the same density and velocity). For model 02, the best-fitting isothermal equivalent has 
a kinetic temperature of 11200K (about 60% of T^s), and its comparison to the radiative 
equilibrium model is shown in Figure 10. The middle panel, which shows the ratio between 
the radiative equilibrium and isothermal IR excesses, indicates how well the isothermal 
model reproduces the full radiative equilibrium solution. Surprisingly, there is little deviation 
(< 20%) for all inclination angles. In detail, the excesses are nearly identical for both 
short and long wavelengths, while the largest differences occur for intermediate wavelengths. 



The reason an isothermal model works so well at very long wavelengths (A > 100 /xm) 
is because the disk is so optically thick that its pseudo-photosphere occurs in the isothermal 
region. For the pole-on case, the effective radius, it!efr, of the pseudo-photosphere is the 
location where the vertical optical depth TA(-Reff) = 1- It is straightforward to show (see 
Appendix A, eq. [A15]) that 
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1 < A < 100 //m. 



i?efr(A) oc A2"-/3-3V2, 



(15) 
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where s and u are the exponents of the power-laws describing the temperature and Gaunt 
factors, respectively (see eqs. [Al] and [A2]). Using /3 = 1.5, n — 3.5, s — 0, and u — 0.23, 
we find (for model 02) that 

i?eff(A) = 1.27R,{X/lijmf^\ (16) 

For A = 1000 /um, this gives Rcs ~ 21i?^, verifying that the effective photospheric radius is 
indeed in the isothermal region for large wavelengths. 

At shorter wavelengths (A < 100 //m), the effective radius is in the non- isothermal re- 
gion. As a result, the isothermal IR excess begins to differ from the radiative equilibrium 
excess. These differences directly reflect the temperature structure of the disk midplane. Us- 
ing eq. (16) to translate from R^s to wavelength, we plot (in the bottom panel of Fig. 10) the 
midplane temperature T(i?cfr) as a function of wavelength. Note the similarity between the 
flux departures (middle panel) and the midplane temperatures (bottom panel). In particu- 
lar, the largest difference occurs near the temperature minimum of the disk. This suggests 
the interesting, albeit technically challenging, possibility of inverting the observed IR excess 
to obtain the midplane temperature profile. The flux excess in the optically thick limit (c.f., 
eq. [A 17]), is given approximately by 

Zx-l ^ 7lBxiRes)Rls/Fx : 
^ T{Res) f -Reff V 

Consequently, the excess ratio 

{Zx - 1)^^ n^f^ 

(Z, - 1)-° " Ti,o V K 

which demonstrates why the flux departures track the the midplane temperatures. Although 
it appears that the flux departure depends linearly on the midplane temperature, one must 
recall that i?efr is temperature-dependent via the opacity. As the midplane radiative equi- 
librium temperature decreases, the opacity increases (see eq. [A5]), making R^ larger than 
the corresponding isothermal value. This increase in the effective area of the disk pseudo- 
photosphere partially offsets the decrease of disk source function, Bx{Res)- The net result 
is that the flux departure is smaller than naively expected. Thus detailed modeling will be 
required to invert the flux departures to obtain the disk midplane temperature versus radius. 

The above discussion explains the flux differences in the 1-100 iim range. It does not, 
however, explain why the isothermal and radiative equilibrium models agree for wavelengths 
A < 1 /im, where the midplane temperature continues to rise well above the isothermal model 
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(see bottom panel of Fig. 10). To help understand this result, we show the wavelength- 
behavior of the IR excess versus temperature in Figure 11. As discussed previously, the IR 
excess increases with temperature for optically thick wavelengths. However for optically thin 
wavelengths, the excess decreases owing to the smaller f-f cmissivity at lower temperatures 
(see eq. [A13]). As a consequence, marginally thick wavelengths have little temperature 
dependence, as shown by the 1 /xm curve (dotted line) in Figure 11. Since {Zx — 1) is a weak 
function of T in the 0.5-1 //m range, isothermal models reproduce the flux excess quite well 
in this region. 

In general, wc have seen that radiative equilibrium models are well-approximated by an 
equivalent isothermal model with a temperature Tigo = 0.6Tefj. However, the temperature 
structure is not the only factor influencing the IR SED. As shown in the appendix, the 
slope of the excess depends on the disk flaring parameter /3 (see eq. [A17]). Furthermore, 
as may be seen in Figures 9 and 10, the slope of the IR excess changes with inclination 
angle (and wavelength). The reason for the inclination-dependence is that, as i?efr changes 
with wavelength, the shape of the flared disk pseudo-photosphere is not self-similar (except 
for the pole-on case). In contrast for an unflared disk {(3 = 1), the shape is self-similar 
(with wavelength) for all inclination angles, so the slope of the IR excess is independent 
of inclination, as demonstrated in Figure 12. We conclude that the primary limitation of 
Waters' method for determining the disk density distribution from the slope of the IR excess 
is not the isothermal assumption; it is the assumption of a particular geometry (e.g., a 
constant opening angle) for the disk. 



6. Discussion and Conclusion 

In this paper, we present a new 3-D NTLE Monte Carlo code that fully solves the radi- 
ation transfer and radiative equilibrium for arbitrary gas density and velocity distributions. 
Even though our method was developed independently, it is similar in many respects to the 
Monte Carlo transition probabihty method of Lucy (2002, 2003). However, our method dif- 
fers from Lucy's in a fundamental way, which is that our photon absorption and reemission 
mechanisms are uncorrelated — we sample the entire emissivity distribution regardless of 
the absorption process. This allows us to dispense with Lucy's macro atoms, along with 
their associated internal transitions and Monte Carlo transition probabilities. Tests of our 
code against an independent (1-D spherically symmetric) NLTE code show that our method 
rapidly converges to the correct results, typically requiring only a few iterations. This rapid 
convergence also occurs for the non-trivial 2-D disk geometry investigated here. 

Since the primary purpose of our NTLE Monte Carlo method is the study of gaseous 
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circumstellar disks, we performed an initial investigation of the temperature and ionization 
structure of Be star disks. Our results for the disk temperature show that its behavior is 
a hybrid between that of optically thick disks of Young Stellar Objects and optically thin 
winds from Hot Stars. We find that the temperature initially drops very quickly close to 
the stellar photosphere, and the temperature profile is well-described by a fiat blackbody 
reprocessing disk (Adams et al. 1987). When the disk becomes optically thin vertically, the 
temperature departs from this curve and rises back to the optically thin radiative equilibrium 
temperature, which is approximately constant in the winds of Hot Stars. We find that the 
average temperature in these isothermal regions is about 60% of Teflf, which is cooler than 
the 80% value frequently quoted in the literature (e.g.. Waters 1986). Similarly, the upper 
layers of the disk heat up to the optically thin radiative equilibrium temperature, which 
again is approximately isothermal. This behavior is analogous to the hot upper layers of 
YSO disks (see Chiang & Goldrcich 1997). However, unlike YSO disks whose flaring controls 
their heating, we find that the disk temperature of Be star disks is infiuenced very little by 
the disk flaring. 

As to be expected, we flnd that the disk is entirely ionized, but there is a region close 
to the midplane (from about 3-6 R^, outwards) for which the neutral hydrogen fraction 
ranges from a few % up to 50%, depending on the model density. This site of low ionization 
coincides with the temperature minimum in the disk midplane. Because the location of 
the temperature minimum is controlled by the vertical optical depth, we find that as the 
density increases, the temperature minimum occurs at larger radii, and the disk becomes 
more neutral owing to the decrease in photoionization at high optical depths. This suggests 
the temperature minimum also is likely to be the site for the formation of low ionization 
state lines (such as Fe II) that produce the shell features observed in Be stars. 

Finally, to investigate the consequences of the non-isothermal temperature structure on 
the IR excess of Be stars, we calculated the SEDs for our models and compared them to 
equivalent isothermal models. We find the effects of the temperature structure are small 
(< 20% change in the flux), and that the IR excess is well-described by an isothermal model 
whose temperature is about 0.6Tcfj. Hence the slope of the IR excess is less sensitive to the 
disk temperature structure than it is to the disk geometry (e.g., inclination angle and flaring). 
We conclude the main hmitation for using the IR excess to determine disk properties is the 
use of an assumed disk geometry. 

In the next paper of this series, we will address the question of disk geometry by self- 
consistently determining both the disk temperature and density. To do so, we will use the 
radiative equilibrium temperatures to solve the vertical hydrostatic equilibrium structure 
for a Keplerian disk with steady-state viscous accretion (or decretion). This self-consistent 
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solution for the disk structure and radiative transfer will permit us to quantitatively test 
whether or not the SEDs and intrinsic polarization of Be stars can be produced by a viscous 
decretion disk. 
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A. IR Excess of Pole-on Flared Disks 

Following a method similar to that employed by Waters (1986), we consider a star 
surrounded by a fully ionized disk with density given by eqs. (9) and (10). Let us assume a 
power-law distribution for the temperature 

T{w) = ToW-' , (Al) 

where w — zu/R^,, as well as a power- law for the free-free plus bound-free Gaunt factors 

9. = *(^)". (A2) 

Since we assume the disk is viewed pole-on, the intensity for < 1 is not changed by 
the presence of the disk. Hence, the intensity may be written as 

r(^^-S ^A(^eff) (^<1), 

^'^^^ ~ \ BfiT) [1 - e--^(^)] (tu > 1), ^^^^ 

where we assumed, for simplicity, that the disk is vertically isothermal and that the star 
emits a blackbody spectrum (for the spectral range of interest, A > 1 /xm, a hot star emits 
a spectrum that is similar to a blackbody). 

The optical depth along the line of sight is 

/oo 
Kx{w)dz, (A4) 
■oo 



where the opacity 



Kx = 3.69 lO^T-V^c-^A^ (1 - e-^'/^^^) gxX^- (A5) 
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For long wavelengths and the typical temperature range of Be star disks (several thousand 
K) , the opacity can be written as 

= Kot^''/'"'"e-^'/^'A2+", (A6) 

with 

Ko = 1.77 lO-^To-^/'c-^^oAo"^. (A7) 
Performing the integration in eq. (A4), we find 

T,{w) = r;^3^/2-2"+^ (A8) 

with t'^ = 0Fko^oA^+". 

The total flux (star plus disk) is given by 

Fx = / d(f) Ix{'^)wdw. (A9) 
Jo Jo 

Defining the excess flux Z\ — \ = {Fx — F^)/F^, we have 

- 1 = 2 || (l - ^ du^ , (AlO) 

which can be written in the Rayleigh- Jeans hmit as 

Z,-l^2^ r ( 1 - e-^A^"^^"''-^") d^. (All) 

TeS Jl ^ ^ 

Equation (All) can be integrated numerically to obtain the fiux excess. However in the 
limiting cases of optically thin and optically thick disks, useful analytical approximations 
can be obtained. In the optically thin case, 1 — exp (— Tj^ro^*/^~^"'+^) fa Tj^ro^*/^~^"+^, so 
eq. (All) becomes 

2^-^, -t' (A12) 

^ reff2 + s/2 + /?-2n ^ ^ ' 

Because r'^ oc T^^^^, it follows that the IR excess in the optically thin case will scale with To 



as 



Zx-l<x T-^'\ (A13) 



In the optically thick limit, a useful approximation is to consider that most of the 
radiation comes from a pseudo-photosphere for which the effective radius ^eff(A) is defined 
so that 

TA(^efr) = 1- (A14) 
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Using eq. ( A8) , we determine that 



u+2 

_ ' \ \ 2n-/3-3s/2 

-Refr(A) = -Ro ( — 



(A15) 



where 



Rq — 



R. 



With i?efr so defined, the fiux excess for the opticaUy thick disk is 



(A16) 



J-eS J I 



2 Tn 



2-sT, 



eflf 



2-s 



(tt+2)(2-«) 
\ \ 2n-/3-3s/2 



(A17) 
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interml trmtsMons 

Fig. 1. — Photon-Macro Atom Interaction. In Lucy's method (Lucy 2002) when a photon 
(r-packet) is absorbed, it either activates the macro atom {A A*), or it is converted into 
a /c-packet (free electron energy). If it activates the macro atom, the macro atom undergoes 
internal transitions (possibly producing a /c-packet and deactivating) until a photon packet is 
eventually reemitted {A* — > A). If a A;-packet is produced, the fc-packet can either reactivate 
the macro atom, or it can emit a photon packet directly. In any case, the absorption of a 
photon eventually leads to the reemission of a new photon. 
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Fig. 2. — Stability Test. The electron temperature (level corresponds to Tg/lO^") and 
hydrogen level populations (up to level 10) are shown for 30 successive iterations. The left 
panel shows the results (triangles) for N — 10^ PPs, while the right panel is for N = 10^ 
PPs. Note that the variation of the values is centered on the correct (LTE) values (stars). 
The R.M.S. variation, e, for various levels is indicated. 
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Fig. 3. — Convergence Test. The hydrogen ground level population is shown for a sequence 
of 30 iterations (Diamonds: N — 10^; Squares: N — 10^; Circles: N — 10^). Iteration zero 
refers to the initial population. The straight line shows the correct level population. 
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Fig. 4. — NLTE Validation. Comparison between the results of the MacFarlane et al. 
(1993) code (open circles) and our code (filled circles). Shown are the n = 2 hydrogen level 
populations as a function of distance from the star. From bottom to top, the curves are for 
logM = —8.7, —7.7 and —6.7, respectively. 
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Fig. 5. — Hydrogen Level Populations. Shown are the level populations for n = 1 (black), 
n = 2 (blue), and n = 3 (green) as a function of radius. The solid curves are for model 01 
(low density) while the dotted curves are for model 04 (high density). Each panel shows the 
populations at a given latitude (angle above midplane), as indicated. 
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Fig. 6. — Disk Temperature. Shown is the temperature structure for model 02 at three radial 
scales. The dashed lines correspond to the curves \z\ = H{zu). 
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Fig. 7. — Mid-plane Temperature. The radial dependence of the equatorial temperature is 
shown for models 01, 02, 03, and 04 (increasing density, respectively). The thick curve is a 
fit, for vj between and ^dep, to a flat blackbody reprocessing disk (Adams et al. 1987). 
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Fig. 8. — Effects of Disk Flaring. The temperature versus radius is shown for six values 
of the disk flaring exponent, (3 = 1.0-1.5 in steps of 0.1 (black, blue, green, red, orange, 
and purple lines, respectively). The different panels show the temperature along rays with 
a given latitudelatitude (left) and for a given fraction of the disk scaleheight (right). 
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Fig. 9. — Spectral Energy Distribution. Shown is the emergent flux for Model 02, viewed at 
three different inclinations, as indicated. The dash-dotted line indicates the input flux from 
the star (a 19000 K Kurucz model atmosphere). 
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Fig. 10. — IR Excess for Model 02. The top panel compares the radiative equilibrium 
IR excess (solid lines) to that of an equivalent isothermal model (dashed lines) for three 
inclinations: pole-on (black), 60° (red), and edge-on (blue). The isothermal model has a 
temperature Ti^o = O.GTes, which yields the best fit. The middle panel shows the ratio of 
the IR excesses for the radiative equilibrium and isothermal models, while the bottom panel 
shows the disk mid-plane temperature at i?efr(A), the radius of the pseudo-photosphere as a 
function of wavelength (see text). Note that the deviation between the radiative equilibrium 
and isothermal excesses roughly tracks the shape of the mid-plane temperature. 
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Fig. 11. — Temperature-dependence of the IR Excess. The IR excess as a function of 
temperature is shown for different wavelengths, as indicated. Each curve was normalized 
to the excess at 5000 K. Note that, for the optically thin (shorter) wavelengths, the excess 
decreases with increasing temperature, while for the optically thick (longer) wavelengths, 
the excess increases. 
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Fig. 12. — Infrared excess vs. wavelength for model 05 {(3 = 1.0). Shown are results for 
three inclinations: pole-on (black), 60° (red), and edge-on (blue). 
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Table 1. Adopted Stellar Parameters. 



Parameter Value 



Spectral Type 
R. 

^eff 

logs' 



t'crit 



Table 2. Disk Model Parameters. 



Model 


Po 














(g cm~'^ 






(K) 






01 


1.66 10" 


-11 


1.5 


17600 


1.7 


.063 


02 


4.15 10- 


-11 


1.5 


19100 


2.0 


.081 


03 


8.39 10- 


-11 


1.5 


19600 


2.5 


.072 


04 


1.66 10- 


-10 


1.5 


20500 


3.0 


.081 


05 


4.15 10- 


-11 


1.0 


18600 


2.0 


.083 


06 


4.15 10- 


-11 


1.1 


18800 


2.0 


.078 


07 


4.15 10" 


-11 


1.2 


18200 


2.0 


.077 


08 


4.15 10" 


-11 


1.3 


18700 


2.0 


.080 


09 


4.15 10" 


-11 


1.4 


18600 


2.0 


.085 



B3 IV 

5.6Rq 
19000 K 
4 

9.2Mo 
560 km s~^ 



